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Abstract 

We study the problem of Robust Least Squares Regression (RLSR) where several response variables 
can be adversarially corrupted. More specifically, for a data matrix X £ and an underlying model 

w*, the response vector is generated as y = X^w* + b where b £ R" is the corruption vector supported 
over at most C -n coordinates. Existing exact recovery results for RLSR focus solely on Li-penalty based 
convex formulations and impose relatively strict model assumptions such as requiring the corruptions b 
to be selected independently of X. 

In this work, we study a simple hard-thresholding algorithm called Torrent which, under mild 
conditions on X, can recover w* exactly even if b corrupts the response variables in an adversarial 
manner, i.e. both the support and entries of b are selected adversarially after observing X and w*. Our 
results hold under deterministic assumptions which are satisfied if X is sampled from any sub-Gaussian 
distribution. Finally unlike existing results that apply only to a fixed w*, generated independently of X, 
our results are universal and hold for any w* £ R^. 

Next, we propose gradient descent-based extensions of Torrent that can scale efficiently to large 
scale problems, such as high dimensional sparse recovery and prove similar recovery guarantees for these 
extensions. Empirically we Hnd Torrent, and more so its extensions, offering significantly faster recovery 
than the state-of-the-art Li solvers. For instance, even on moderate-sized datasets (with p = 50K) with 
around 40% corrupted responses, a variant of our proposed method called Torrent-HYB is more than 
20 X faster than the best Li solver. 


“If among these errors are some which appear too large to be admissible, 
then those equations which produced these errors will be rejected, as com¬ 
ing from too faulty experiments, and the unknowns will be determined by 
means of the other equations, which will then give much smaller errors. ” 
A. M. Legendre, On the Method of Least Squares. 1805. 


1 Introduction 


Robust Least Squares Regression (RLSR) addresses the problem of learning a reliable set of regression 
coefficients in the presence of several arbitrary corruptions in the response vector. Owing to the wide- 
applicability of regression, RLSR features as a critical component of several important real-world applications 
in a variety of domains such as signal processing |13j , economics [12] , computer vision [161115] , and astronomy 

[Hi- 

Given a data matrix X = [xi,... ,x„] with n data points in K*’ and the corresponding response vector 
y £ R", the goal of RLSR is to learn a w such that. 
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That is, we wish to simultaneously determine the set of corruption free points S and also estimate the 
best model parameters over the set of clean points. However, the optimization problem given above is non- 
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convex (jointly in w and S) in general and might not directly admit efficient solutions. Indeed there exist 
reformulations of this problem that are known to be NP-hard to optimize m- 

To address this problem, most existing methods with provable guarantees assume that the observations 
are obtained from some generative model. A commonly adopted model is the following 


y = x'^w* + b, (2) 

where w* G is the true model vector that we wish to estimate and b G K” is the corruption vector that 
can have arbitrary values. A common assumption is that the corruption vector is sparsely supported i.e. 
|lb||o < a • n for some a > 0. 

Recently, m and [9] obtained a surprising result which shows that one can recover w* exactly even when 
a 1, i.e., when almost all the points are corrupted, by solving an Li-penalty based convex optimization 
problem: minw,b ||w||j+A ||b||^, s.t., A^w+b = y. However, these results require the corruption vector b to 
be selected oblivious of X and w*. Moreover, the results impose severe restrictions on the data distribution, 
requiring that the data be either sampled from an isotropic Gaussian ensemble |15j . or row-sampled from 
an incoherent orthogonal matrix [9]. Finally, these results hold only for a fixed w* and are not universal in 
general. 

In contrast, [3] studied RLSR with less stringent assumptions, allowing arbitrary corruptions in response 
variables as well as in the data matrix A, and proposed a trimmed inner product based algorithm for the 
problem. However, their recovery guarantees are significantly weaker. Firstly, they are able to recover w* 
only upto an additive error a^/p (or a^/s if w* is s-sparse). Hence, they require a < y/p just to claim a 
non-trivial bound. Note that this amounts to being able to tolerate only a vanishing fraction of corruptions. 
More importantly, even with n —> oo and extremely small a they are unable to guarantee exact recovery 
of w*. A similar result was obtained by [S], albeit using a sub-sampling based algorithm with stronger 
assumptions on b. 

In this paper, we focus on a simple and natural thresholding based algorithm for RLSR. At a high level, 
at each step t, our algorithm alternately estimates an active set St of “clean” points and then updates the 
model to obtain by minimizing the least squares error on the active set. This intuitive algorithm seems 
to embody a long standing heuristic first proposed by Legendre [T] over two centuries ago (see introductory 
quotation in this paper) that has been adopted in later literature [TUI El] as well. However, to the best of 
our knowledge, this technique has never been rigorously analyzed before in non-asymptotic settings, despite 
its appealing simplicity. 

Our Contributions: The main contribution of this paper is an exact recovery guarantee for the thresh¬ 
olding algorithm mentioned above that we refer to as Torrent-FC (see Algorithm [^. We provide our 
guarantees in the model given in where the corruptions b are selected adversarially but restricted to have 
at most a ■ n non-zero entries where a < 1/2 is a global constant dependent only on aQ Under deter¬ 
ministic conditions on A, namely the subset strong convexity (SSC) and smoothness (SSS) properties (see 
Definition [^ , we guarantee that Torrent-FC converges at a geometric rate and recovers w* exactly. We 
further show that these properties (SSC and SSS) are satisfied w.h.p. if a) the data A is sampled from a 
sub-Gaussian distribution and, b) n > p\ogp. 

We would like to stress three key advantages of our result over the results of (TU] [U] : a) we allow b to be 
adversarial, i.e., both support and values of b to be selected adversarially based on A and w*, b) we make 
assumptions on data that are natural, as well as significantly less restrictive than what existing methods 
make, and c) our analysis admits universal guarantees, i.e., holds for any w*. 

We would also like to stress that while hard-thresholding based methods have been studied rigorously for 
the sparse-recovery problem Hi, hard-thresholding has not been studied formally for the robust regression 
problem. Moreover, the two problems are completely different and hence techniques from sparse-recovery 
analysis do not extend to robust regression. 

^Note that for an adaptive adversary, as is the case in our work, recovery cannot be guaranteed for a > 1/2 since the 
adversary can introduce corruptions as = xT (w — w*) for an adversarially chosen model w. This would make it impossible 
for any algorithm to distinguish between w* and w thus making recovery impossible. 
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Despite its simplicity, Torrent-FC does not scale very well to datasets with large p as it solves least 
squares problems at each iteration. We address this issue by designing a gradient descent based algorithm 
(Torrent-GD), and a hybrid algorithm (TORRENX-Hyb), both of which enjoy a geometric rate of con¬ 
vergence and can recover w* under the model assumptions mentioned above. We also propose extensions 
of Torrent for the RLSR problem in the sparse regression setting where p ^ n but ||w*||g = s* <C p. 
Our algorithm Torrent-HD is based on Torrent-FC but uses the Iterative Hard Thresholding (IHT) 
algorithm, a popular algorithm for sparse regression. As before, we show that Torrent-HD also converges 
geometrically to w* if a) the corruption index a is less than some constant C, b) X is sampled from a 
sub-Gaussian distribution and, c) n > s* logp. 

Finally, we experimentally evaluate existing Li-based algorithms and our hard thresholding-based algo¬ 
rithms. The results demonstrate that our proposed algorithms (Torrent-(FG/GD/HYB)) can be signifi¬ 
cantly faster than the best Li solvers, exhibit better recovery properties, as well as be more robust to dense 
white noise. For instance, on a problem with 50K dimensions and 40% corruption. Torrent- HYB was 
found to be 20 x faster than Li solvers, as well as achieve lower error rates. 

Paper Organization: We give a formal definition of the RLSR problem in the next section. We then 
introduce our family of algorithms in Section and prove their convergence guarantees in Section We 
present extensions to sparse robust regression in Section]^ and empirical results in Section]^ 


2 Problem Formulation 

Given a set of data points X = [xi,X 2 ,... ,x„], where x^ G and the corresponding response vector 
y S K", the goal is to recover a parameter vector w* which solves the RLSR problem Q. We assume that 
the response vector y is generated using the following model: 

y = y* -I- b -I- e, where y* = A^w*. 

Hence, in the above model, Q reduces to estimating w*. We allow the model w* representing the regressor, 
to be chosen in an adaptive manner after the data features have been generated. 

The above model allows two kinds of perturbations to yi - dense but bounded noise Si (e.g. white noise 
Ei ~ A/'(0, (T^), CT > 0), as well as potentially unbounded corruptions ft - to be introduced by an adversary. 
The only requirement we enforce is that the gross corruptions be sparse. 

e shall represent the dense noise vector, for example e ~ A/'(0, • Inxn), and b, the corruption vector 

such that Ijbllp < a • n for some corruption index a > 0. We shall use the notation S'* = supp(b) C [n] 
to denote the set of “clean” points, i.e. points that have not faced unbounded corruptions. We consider 
adaptive adversaries that are able to view the generated data points x^, as well as the clean responses y* 
and dense noise values Ei before deciding which locations to corrupt and by what amount. 

We denote the unit sphere inp dimensions using S^~^. For any 7 G (0,1], we let Sj = {S C [n] : |S| = 7 • n} 
denote the set of all subsets of size 7 • n. For any set S, we let Xs := G denote the matrix 

whose columns are composed of points in that set. Also, for any vector v G M" we use the notation vg 
to denote the |S'|-dimensional vector consisting of those components that are in S. We use Ai„in(A) and 
Amax(A) to denote, respectively, the smallest and largest eigenvalues of a square symmetric matrix X. We 
now introduce two properties, namely. Subset Strong Convexity and Subset Strong Smoothness, which are 
key to our analyses. 

Definition 1 (SSC and SSS Properties). A matrix X G satisfies the Subset Strong Convexity Property 
(resp. Subset Strong Smoothness Property^ at level 7 with strong convexity constant (resp. strong 
smoothness constant A^) if the following holds: 

X-y min Aniin(AgAg ) ^ mt^Amax(AgAg ) ^ A.y. 

Remark 1. We note that the uniformity enforced in the definitions of the SSC and SSS properties is not for 
the sake of convenience but rather a necessity. Indeed, a uniform bound is required in face of an adversary 
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Algorithm 1 Torrent: Thresholding Operator-based Algorithm 3 Torrent-GD 
Robust RegrEssioN meThod 


Input: Training data {xi, yi} ,i = 1... n, step length tj, 
thresholding parameter /3, tolerance e 


w° ^ 0, S'o = [n], t ^ 0, r° ■<- y 
while Iks II > e do 


Input: Current model w, current active set S, step size rj 
1: g t- As(AJw - ys) 

2: return w — r] ■ g 


w*+i 
r‘+i 4 
St+i ^ HT(r 


UPDATE(w*,St,77,r‘,5t-i) 
{Vi - (w*+kxi)) 


Algorithm 4 Torrent-HYB 


t+i 




Input: Current model w, current active set S, step size 
rj, current residuals r, previous active set S' 

1: // Use the GD update if the active set S is 


7: end while 


changing a lot 

8: return w* 

2 

if IS’VS"! > A then 


3 

w' ^ UPDATE-GD(w, S, g, r, S') 


- 4 

else 

Algorithm 2 Torrent-FC 

5 

//If stable, use the FC update 

Input: Current model w, current active set S 

6 

w' c- UPDATE-FC(w, S) 

1: return arg min {yi — {w, Xi))^ 

7 

end if 

iGS 

8 

return w' 


which can perforin corruptions after data and response variables have been generated, and choose to corrupt 
precisely that set of points where the SSC and SSS parameters are the worst. 

3 Torrent: Thresholding Operator-based Robust Regression Method 

We now present Torrent, a Thresholding Operator-based Robust RegrEssioN meThod for performing 
robust regression at scale. Key to our algorithms is the Hard Thresholding Operator which we define below. 

Definition 2 (Hard Thresholding Operator). For any vector v S M", let G Sn he the permutation that 
orders elements of v in ascending order of their magnitudes i.e. |v£,.^(i)| < |vct^( 2 )| < ... < |vo.^(„)|. Then 
for any k <n, we define the hard thresholding operator as 

HT(v;A:) = {i G [n] : cr~^{i) < k} 

Using this operator, we present our algorithm Torrent (Algorithm]^ for robust regression. Torrent 
follows a most natural iterative strategy of, alternately, estimating an active set of points which have the 
least residual error on the current regressor, and then updating the regressor to provide a better fit on this 
active set. We offer three variants of our algorithm, based on how aggressively the algorithm tries to fit the 
regressor to the current active set. 

We first propose a fully corrective algorithm Torrent-FC (Algorithm]^ that performs a fully corrective 
least squares regression step in an effort to minimize the regression error on the active set. This algorithm 
makes significant progress in each step, but at a cost of more expensive updates. To address this, we then 
propose a milder, gradient descent-based variant Torrent-GD (Algorithm]^ that performs a much cheaper 
update of taking a single step in the direction of the gradient of the objective function on the active set. 
This reduces the regression error on the active set but does not minimize it. This turns out to be beneficial 
in situations where dense noise is present along with sparse corruptions since it prevents the algorithm from 
overfitting to the current active set. 

Both the algorithms proposed above have their pros and cons - the FC algorithm provides significant 
improvements with each step, but is expensive to execute whereas the GD variant, although efficient in exe¬ 
cuting each step, offers slower progress. To get the best of both these algorithms, we propose a third, hybrid 
variant Torrent- HYB (Algorithm]^ that adaptively selects either the FC or the GD update depending on 
whether the active set is stable across iterations or not. 
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In the next section we show that this hard thresholding-based strategy offers a linear convergence rate 
for the algorithm in all its three variations. We shall also demonstrate the applicability of this technique to 
high dimensional sparse recovery settings in a subsequent section. 


4 Convergence Guarantees 


For the sake of ease of exposition, we will first present our convergence analyses for cases where dense noise 
is not present i.e. y = -|- b and will handle cases with dense noise and sparse corruptions later. 

We first analyze the fully corrective Torrent-FC algorithm. The convergence proof in this case relies on 
the optimality of the two steps carried out by the algorithm, the fully corrective step that selects the best 
regressor on the active set, and the hard thresholding step that discovers a new active set by selecting points 
with the least residual error on the current regressor. 

Theorem 3. Let X = [xi,... ,x„] G be the given data matrix and y = -|- b &e the corrupted 

output with |jb||o < a ■ n. Let Algorithm^ be executed on this data with the thresholding parameter set to 
(3 > a. Let Sq invertible matrix such that X = satisfies the SSC and SSS properties at level 

. If the data satisfies < 1; then after 

t = O ^log j j iterations, Algorithm^ obtains an e-accurate solution w* i.e. ||w* — w*||2 < e. 

Proof (Sketch). Let r‘ = y — X^w* be the vector of residuals at time t and Ct = Xs^Xj^. Also let 
S'* = supp(b) be the set of uncorrupted points. The fully corrective step ensures that 

= Cfi^Xs.ys, = Cfi^Xs, (Xjw* + bsj = w* + Cfi^Xs.hs,, 


7 with constants and A.y respectively (see Definition 
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whereas the hard thresholding step ensures that 


..i+l 


< 




- Ids. 


2 

2 ' 


Combining the two gives us 


PS't+i||2 < 


X 


s.\St+iC't ^ + 2 • ^Xsjbs, 

2 


Cl 








\St+i ^Stbst ^ + (^Xs^Xj^) ^Stbs^ 

bsJIs llbsj+J^ ’ 


|b5jl2 + 2 


Ai 


where Ci follows from setting X = Sq ^^^X and X^Cfi^Xs' = Xj(XsjXjJ“^Xs' and (2 follows from the 
SSC and SSS properties, IjbsJlQ < ||b||p < (3 ■ n and |S*\St+i| < (3 ■ n. Solving the quadratic equation and 
performing other manipulations gives us the claimed result. □ 

Theorem relies on a deterministic {fixed design) assumption, specifically < 1 in order to 

guarantee convergence. We can show that a large class of random designs, including Gaussian and sub- 
Gaussian designs actually satisfy this requirement. That is to say, data generated from these distributions 
satisfy the SSC and SSS conditions such that ^ ^ probability. Theorem explicates this 

for the class of Gaussian designs. 

Theorem 4. Let X = [xi,...,x„] G be the given data matrix with each x^ ~ A/'(0, S). Let y = 

X^w* -t- b and ||b||o < a ■ n. Also, let a < (3 < and n > Q (p-|-logy). Then, with probability at 

least 1 — i5, the data satisfies ^ TO’ More specifically, after T > 10log j iterations of 

Algorithm^with the thresholding parameter set to (3, we have ||w^ — ■w*|| < e. 
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Remark 2. Note that Theorem 4 provides rates that are independent of the condition number of the 

I_I 

distribution. We also note that results similar to Theorem|^can be proven for the larger class of sub-Gaussian 
distributions. We refer the reader to Section for the same. 

Remark 3. We remind the reader that our analyses can readily accommodate dense noise in addition to 
sparse unbounded corruptions. We direct the reader to Appendix which presents convergence proofs for 
our algorithms in these settings. 

Remark 4. We would like to point out that the design requirements made by our analyses are very mild when 
compared to existing literature. Indeed, the work of [15] assumes the Bouquet Model where distributions 
are restricted to be isotropic Gaussians whereas the work of [S] assumes a more stringent model of sub¬ 
orthonormal matrices, something that even Gaussian designs do not satisfy. Our analyses, on the other 
hand, hold for the general class of sub-Gaussian distributions. 

We now analyze the Torrent-GD algorithm which performs cheaper, gradient-style updates on the 
active set. We will show that this method nevertheless enjoys a linear rate of convergence. 

Theorem 5. Let the data settings be as stated in Theorem^ and let Algorithm^ be executed on this data 
with the thresholding parameter set to P > a and the step length set to rj = -. If the data satisfies 

max 1 — ryAi-^g} < then after t = O ^log e)) iterations, obtains an e-accurate 

solution w‘ i.e. ||w* — w *||2 < e. 

Similar to Torrent-FG, the assumptions made by the Torrent-GD algorithm are also satisfied by the 
class of sub-Gaussian distributions. The proof of Theorem given in Appendix [D| details these arguments. 
Given the convergence analyses for Torrent-FG and GD, we now move on to provide a convergence analysis 
for the hybrid Torrent- HYB algorithm which interleaves FG and GD steps. Since the exact interleaving 
adopted by the algorithm depends on the data, and not known in advance, this poses a problem. We address 
this problem by giving below a uniform convergence guarantee, one that applies to every interleaving of the 
FG and GD update steps. 

Theorem 6. Suppose Algorithm^is executed on data that allows Algorithms^an a convergence rate of 
ijFC ci'iT'd rjGD respectively. Suppose we have 2 • rjpc ' Vgb < 1- Then for any interleavings of the FC and GD 
steps that the policy may enforce, after t = O ^log iterations, AZgorzt/im ensures an e-optimal 

solution i.e. ||w‘ — w*|| < e. 

We point out to the reader that the assumption made by Theorem Hi.e. 2 • Pfo ■ ??GD < 1 is readily 
satisfied by random sub-Gaussian designs, albeit at the cost of reducing the noise tolerance limit. As we 
shall see, Torrent-HYB offers attractive convergence properties, merging the fast convergence rates of the 
FG step, as well as the speed and protection against overhtting provided by the GD step. 

5 High-dimensional Robust Regression 

In this section, we extend our approach to the robust high-dimensional sparse recovery setting. As before, 
we assume that the response vector y is obtained as: y = Ai^w* + b, where |jb||o < a ■ n. Bowever, this 
time, we also assume that w* is s*-sparse i.e. ||w*||q < s*. 

As before, we shall neglect white/dense noise for the sake of simplicity. We reiterate that it is not possible 
to use existing results from sparse recovery (such as [3l|6]) directly to solve this problem. 

Our objective would be to recover a sparse model w so that ||w — w *||2 < e. The challenge here is to 
forgo a sample complexity of n > p and instead, perform recovery with n ~ s* logp samples alone. For this 
setting, we modify the FG update step of Torrent-FG method to the following: 

^ II ~ ’ (3) 
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for some target sparsity level s p. We refer to this modified algorithm as Torrent-HD. Assuming X 
satisfies the RSC/RSS properties (defined below), ([^ can be solved efficiently using results from sparse 
recovery (for example the IHT algorithm [3l[5] analyzed in [ 6 ]). 

Definition 7 (RSC and RSS Properties). A matrix X G will be said to satisfy the Restricted Strong 

Convexity Property (resp. Restricted Strong Smoothness Property^ at level s = S 1 + S 2 with strong convexity 
constant asi+s 2 (i"ssp. strong smoothness constant if the following holds for all ||wi||q < si and 

IIW2II0 < S 2 -' 

Og ||wi - W2II2 < ||A^(wi - W2)||2 < Ls ||wi - W2II2 

For our results, we shall require the subset versions of both these properties. 

Definition 8 (SRSC and SRSS Properties). A matrix X G will be said to satisfy the Subset Restricted 
Strong Convexity (resp. Subset Restricted Strong Smoothness^ Property at level ( 7 , s) with strong convexity 
constant 0 ( 7 ,g) (resp. strong smoothness constant g)^ if for all subsets S G Sj, the matrix Xg satisfies 
the RSC (resp. RSS) property at level s with constant Og (resp. Ls). 


We now state the convergence result for the Torrent-HD algorithm. 

Theorem 9. Let X G be the given data matrix and y = X^w* + b 6 e the corrupted output with 

'lw*||n < s* and ||b||o < a-n. Let Eg be an invertible matrix such that Eg satisfies the SRSC and SRSS 


IQ O UIOUj II U||Q j;;- Cl 10. UC- tOlO OlOVC-l OOUIZ. IIOULI OU. OU-C/t UlOUO A^q 

properties at level ( 7 , 2s + s*) with constants Q;(g., 2 s+s*) ond L(g, 2 g+g*) respectively (see Definition] 


Let 


Algorithm^be executed on this data with the Torrent-HD update, thresholding parameter set top > a, 


and s 


> 32^ 


"(1-/3,2a + , 
“(1-/3,2!,+, 


s) 


If X also satisfies ^ then after t = O ^log iterations, Algorithmj^ obtains 

logp), 


w — w 


l 2 <£- 


e-accurate solution i.e. 

In particular, if X is sampled from a Gaussian distribution Af{0, E) and n> Q ^s* 

for all values of a < (3 < we can guarantee ||w* — w* II2 < e after t = O |^log 
the algorithm (w.p. > 1 — 1/n^^). 


x(S) 

hs) 


then 


iterations of 


Remark 5. The sample complexity required by Theorem]^ is identical to the one required by analyses for 
high dimensional sparse recovery [^, save constants. Also note that Torrent-HD can tolerate the same 
corruption index as Torrent-FC. 


6 Experiments 

Several numerical simulations were carried out on linear regression problems in low-dimensional, as well as 
sparse high-dimensional settings. The experiments show that Torrent not only offers statistically better 
recovery properties as compared to Li-style approaches, but that it can be more than an order of magnitude 
faster as well. 

Data; For the low dimensional setting, the regressor w* G was chosen to be a random unit norm 
vector. Data was sampled as Xi ~ A/’(0, Ip) and response variables were generated as y* = (w*, x^). The set of 
corrupted points S'* was selected as a uniformly random (Q;n)-sized subset of [n] and the corruptions were set 
to bi ~ C (—5 ||y*||g^ , 5 ||y*||g^) for i G S*. The corrupted responses were then generated as Vi = yt + bi +ei 
where "x. A/’(0, cr^). For the sparse high-dimensional setting, supp(w*) was selected to be a random s*-sized 
subset of [p]. Phase-transition diagrams (Figure]^ were generated by repeating each experiment 100 times. 
For all other plots, each experiment was run over 20 random instances of the data and the plots were drawn 
to depict the mean results. 

Algorithms; We compared various variants of our algorithm Torrent to the regularized Li algorithm 
for robust regression mM- Note that the Li problem can be written as min^ ||z|| s.t.Az = y, where 
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Figure 1: (a), (b) and (c) Phase-transition diagrams depicting the recovery properties of the Torrent-FC, Torrent- 
HYB and L\ algorithms. The colors red and bine represent a high and low probability of success resp. A method is 
considered successfnl in an experiment if it recovers w* upto a 10“^ relative error. Both variants of Torrent can 
be seen to recover w* in presence of larger nnmber of corrnptions than the L\ solver, (d) Variation in recovery error 
with the magnitude of corruption. As the corruption is increased, Torrent-FC and Torrent-HYB show improved 
performance while the problem becomes more difficnlt for the Li solver. 

A = [X^ y/mxm] and z* = [w*^ Ab^]^. We used the Dual Augmented Lagrange Multiplier (DALM) Li 
solver implemented by m to solve the Li problem. We ran a fine tuned grid search over the A parameter 
for the Li solver and quoted the best results obtained from the search. In the low-dimensional setting, we 
compared the recovery properties of Torrent-FC (Algorithm]^ and Torrent-HYB (Algorithm]^ with 
the DALM-Li solver, while for the high-dimensional case, we compared Torrent-HD against the DALM-Li 
solver. Both the Li solver, as well as our methods, were implemented in Matlab and were run on a single 
core 2.4GHz machine with 8 GB RAM. 

Choice of Li-solver: An extensive comparative study of various Li minimization algorithms was 
performed by m who showed that the DALM and Homotopy solvers outperform other counterparts both 
in terms of recovery properties, and timings. We extended their study to our observation model and found 
the DALM solver to be signihcantly better than the other Li solvers; see Figure]^ in the appendix. We also 
observed, similar to that the Approximate Message Passing (AMP) solver diverges on our problem as 
the input matrix to the Li solver is a non-Gaussian matrix A = [X'^jl]. 

Evaluation Metric: We measure the performance of various algorithms using the standard L2 error: 
fw = llw —w*|| 2 . For the phase-transition plots (Figure [T), we deemed an algorithm successful on an 
instance if it obtained a model w with error < 10“^ • |lw^| 2 . We also measured the GPU time required 
by each of the methods, so as to compare their scalability. 

6.1 Low Dimensional Results 

Recovery Property: The phase-transition plots presented in Figurej^represent our recovery experiments 
in graphical form. Both the fully-corrective and hybrid variants of Torrent show better recovery properties 
than the Li-minimization approach, indicated by the number of runs in which the algorithm was able to 
correctly recover w* out of a 100 runs. Figurej^shows the variation in recovery error as a function of a in the 
presence of white noise and exhibits the superiority of Torrent-FG and Torrent-HYB over Li-DALM. 
Here again, Torrent-FG and Torrent-HYB achieve significantly lesser recovery error than Li-DALM for 
all a <= 0.5. Figure|^in the appendix show that the variations of ||w — w *||2 with varying p, a and n follow 
a similar trend with Torrent having significantly lower recovery error in comparison to the Li approach. 

Figure m brings out an interesting trend in the recovery property of Torrent. As we increase 
the magnitude of corruption from U (— ||y*||go , ||y*|loo) to U (—20 ||y*||oo , 20 ||y*||oo), the recovery error for 
Torrent-HYB and Torrent-FG decreases as expected since it becomes easier to identify the grossly 
corrupted points. However the Li-solver was unable to exploit this observation and in fact exhibited an 
increase in recovery error. 

Run Time: In order to ascertain the recovery guarantees for Torrent on ill-conditioned problems, we 
performed an experiment where data was sampled as ~ A/’(0, S) where diag(E) ~ 1/(0, 5). Figure plots 

















500 n = 2000 sigma = 0.2 p = 300 n = 1800 alpha = 0.41 kappa = 5 



(a) 


(b) 


p = 10000 n = 2303 s = 50 p = 50000 n = 5410 alpha = 0.4 s = 100 



Figure 2: (a), (b) and (c) Phase-transition diagrams depicting the recovery properties of the Torrent-FC, Torrent- 
HYB and Li algorithms. The colors red and bine represent a high and low probability of success resp. A method is 
considered successfnl in an experiment if it recovers w* npto a 10“^ relative error. Both variants of Torrent can 
be seen to recover w* in presence of larger nnmber of corruptions than the L\ solver, (d) Variation in recovery error 
with the magnitude of corruption. As the corruption is increased, Torrent-FC and Torrent-HYB show improved 
performance while the problem becomes more difficult for the Li solver. 


the recovery error as a function of time. Torrent-HYB was able to correctly recover w* about 50x faster 
than Li-DALM which spent a considerable amount of time pre-processing the data matrix X. Even after 
allowing the Li algorithm to run for 500 iterations, it was unable to reach the desired residual error of 10“"^. 
Figure also shows that our Torrent-HYB algorithm is able to converge to the optimal solution much 
faster than Torrent-FC or Torrent-GD. This is because Torrent-FC solves a least square problem at 
each step and thus, even though it requires significantly fewer iterations to converge, each iteration in itself 
is very expensive. While each iteration of Torrent-GD is cheap, it is still limited by the slow O ((1 — ^)*) 
convergence rate of the gradient descent algorithm, where k is the condition number of the covariance matrix. 
Torrent-BYB, on the other hand, is able to combine the strengths of both the methods to achieve faster 
convergence. 

6.2 High Dimensional Results 

Recovery Property; Figure shows the variation in recovery error in the high-dimensional setting as 
the number of corrupted points was varied. For these experiments, n was set to 5s* log(p) and the fraction of 
corrupted points a was varied from 0.1 to 0.7. While Li-DALM fails to recover w* for a > 0.5, Torrent-BD 
offers perfect recovery even for a values upto 0.7. 

Run Time: Figure shows the variation in recovery error as a function of run time in this setting. 
Li-DALM was found to be an order of magnitude slower than Torrent-BD, making it infeasible for sparse 
high-dimensional settings. One key reason for this is that the Li-DALM solver is significantly slower in 
identifying the set of clean points. For instance, whereas Torrent-BD was able to identify the clean set of 
points in only 5 iterations, it took Li around 250 iterations to do the same. 
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A Convergence Guarantees with Dense Noise and Sparse Corrup¬ 
tions 

We will now present recovery guarantees for the Torrent-FC algorithm when both, dense noise, as well 
as sparse adversarial corruptions are present. Extensions for Torrent-GD and Torrent-HYB will follow 
similarly. 

Theorem 10. Let X = [xi,... ,x„] G be the given data matrix and y = J'f^w* + b + e be the corrupted 

output with sparse corruptions ||b||o < a-n as well as dense bounded noise e. Let Algorithm^be executed on 
this data with the thresholding parameter set to (3 > a. Let Eq be an invertible matrix such that V = Eg 
satisfies the SSC and SSS properties at level 7 with constants Ag, and respectively (see Definition^. If 

the data satisfies < 1? then after t = O ^log iterations, Algorithm 2 obtains an e-accurate 

solution w‘ i.e. ||w* — w *||2 < e + for some constant C > 0. 

Proof. We being by observing that the optimality of the model on the active set St ensures 

||yst = ||xj(w* - w‘+i) +£s, +bsJ 2 < l|y* = lies, +bsj| 2 , 

which, upon the application of the triangle inequality, gives us 


Since ||Xj(w* - w ‘+^)||2 > 

j+i 


(w* — w*' 
' _ I 


^)||2 < 2 ||est +bs. 

we get 


t 112 ■ 


2 2 

|w'- - W- -||2 < - lies, +bsj |2 < 




)• 


The hard thresholding step, on the other hand, guarantees that 

- w*+i) + es,+i + bs,+i 


yst. 


-X 


St+' 


w 


i +1 


< ||ys. - ^S.W*+^||2 

= (w* — + es. 


As before, let CRt+i = S't+i\S'» and MD(+i = St,\St+i. Then we have 


^CRt+i(w* - W + ) + ecR,+i + bcR,+i 


< 


X, 


MD 


t+ 


,(w* - W + ) +eMDt+i 


An application of the triangle inequality and the fact that ||bcR,^.i IL = Il'^S't+i 11 gives us 




^MDt+i(w* - W 


<2v^| |w"' — w 


t+ 1 ' 

t+l|| 


+ 


^CRt+i(w* - W 


t + 1 II 2 
t+1 


^CRt+i 9 ^MDt_Li 


t +1 II 2 


7211, 


l|t>Stll 2 + ( )C~^ + 


<r?.||bs ,||2 + (l + 72)||e|l 


where the second step uses the fact that max{|CRt+i|, iMD^+il} < (d-n and the Cauchy-Schwartz inequality, 
step uses the fac 
12 again gives us 


and the last step uses the fact that for sufficiently small /?, we have 7 := Using the inequality for 

V-^l-ZS 


w — w 


t +11 


12 < 


7 A 1 -/; 


+ I|bstll 2 ) 
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< 


4 + 2^2 

\/^l-/3 


+ 


2-r;‘ 


For large enough n we have > O {\/n), which completes the proof. 

Notice that for random Gaussian noise, this result gives the following convergence guarantee. 


□ 


Corollary 11. Let the date be generated as before with random Gaussian dense noise i.e. y = X'^w*+b + e 
with ||b||o < a-n and e ^ Xf{0, Let Algorithm^be executed on this data with the thresholding parameter 

set to (3 > a. Let Eq be an invertible matrix such that X = XX 

level 7 with constants A.y and respectively (see Definition 


t = 


O ^log iterations, obtains an e-accurate solution 


X satisfies the SSC and SSS properties at 
1). If the data satisfies < 1, then after 

VAi_/3 

12 ^ e + 2tTC, 


where C > Q is the constant in Theorem \l(h 
Proof. Using tail bounds on Chi-squared distributions [7], we get, with probability at least 1 — <5, 


1 


1 


112 < 0 - I n + 2\j n\og - + 2 log - 


Thus, for n > 4log j, we have ||e ||2 < 2an which proves the result. 

Remark 6. We note that the design assumptions made by Theorem 

those made by Theorem]^ and would be satisfied with high probability 
distributions (see Appendixfor details). 


□ 


(i..e < 1) are similar to 

data sampled from sub-Gaussian 


B Proof of Theorem [3] 

Theorem 3. Let X = [xi,... ,x„] G be the given data matrix and y = A^w* -|- b be the corrupted 

output with ||b||o < a • n. Let Algorithmj^be executed on this data with the thresholding parameter set to 
(3 > a. Let Eq be an invertible matrix such that X = Eq satisfies the SSC and SSS properties at level 
7 with constants A.y and A.y respectively (see Definition j^. If the data satisfies < 1, then after 

t = O ^log iterations. Algorithmobtains an e-accurate solution w* i.e. ||w‘ — w *||2 < e. 

Proof. Let r* = y — A^w* be the vector of residuals at time t and Ct = Ag^Aj^. Since Aq, > 0 (something 
which we shall establish later), we get 

w‘+i = Cfi^Xs.yst = w* + hs,) + Cfi^Xs.hs,. 

Thus, for any set S C [n], we have 

=ys - Ajwt+i = bs - AjC("^As^bs^ 


ims, gives 


- Xl^^Ci^Xs,hs,\\l - ||Aj^^^Cr'As,b5,|[ -t- 2 • 




l + 2-hl^^Xl^^Cfi^Xs,hs, 


||b5. - AXCr'As,b5j2 - l^y+.CT'^s.bs, 

^ ||AXCr'As,b5j2 - ||^y+,C'r'^s.b5,|[ + 2 . b^^^A^^^Cr'As.bs, 
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< 

Cs 




Xl 


AS.^^ + 2 • "' ^5.b5. 


C4 -^p 9 -A.P 

<T^-|!b5jl2 + 2-^ 

Ai_/3 Ai-^ 


||bsj|2 llbst^i 


t+i II 2 ' 


where Ci follows since the hard thresholding step ensures 


[.t+i 


— Il’^sl^ll 2 Claim 
^- 1/2 - 


19 


and use the fact 


that P > a), C 2 notices the fact that bs, = 0. ^3 follows from setting X = Eg '^'‘^X and XjC^ ^^ 5 / = 
Xj(XsjXjJ“^Xs/. C 4 follows from the definition of SSC and SSS properties, ||bsj|g < Ijbllg < P ■ n and 
|S'*\S't+i| < P ■ n. Solving the quadratic equation gives us 


bs. 


<(1 + 72).^.libs. 


i + 1 II 2 


Al-/; 


(4) 


Let 77 := denote the convergence rate in Q. We shall show below that for a large family of 

random designs, we have rj < 1 if n > n(p + logj). We now recall from our earlier discussion that 
= w* + Xsjast which gives us 


w*+i-w*|| = ||Cr^Xs.bsJL< 




Ibstlla < 


\/~h 


blU < e, 


- ' Ai_7 

for t > logi ^ ■ Noting that < O establishes the convergence result. 

C Proof of Theorem S] 


□ 


Theorem 4 - Let X = [xi,...,x„] G be the given data matrix with each x^ ~ ^(0, E). Let y = 

X^w* + b and ||b||o < a ■ n. Also, let a < ,0 < A and n > Tl (p + log j). Then, with probability at 

least 1 — A, the data satisfies More specifically, after T > 10 log iterations of 

Algorithmwith the thresholding parameter set to /3, we have ||w^ — w*|| < e. 

Proof. We note that whenever x ^ ^(0, E) then E“^/^x ^ Thus, Theorem 15 assures us that 

with probability at least 1 — (5, the data matrix X = E“^/^X satisfies the SSC and SSS properties with the 
following constants 


A /3 < /3n 1 + 3e,/61og ^ ] + 0 \ \ np + nlog 


A 1-/3 > n - /3n ( 1 + 3eJ61og - ) - fl | y np + nlog - 


Thus, the convergence given be Algorithmj^ when invoked with Eg = E, relies on the quantity p = 

being less than unity. This translates to the requirement (1 + ^/2)Kp < Xi-p. Using the above bounds 
translates that requirement to 


(2 + VT)P ^1 + 36^6 log + o( y I + - log 


(A) 


(S) 
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For n = fl[p + log y), the second quantity (B) can be made as small a constant as necessary. Tackling the 
first quantity [A) turns out to be more challenging. However, we can show that for all /3 < we get 

77 = which establishes the claimed result. Thus, Algorithm can tolerate a corruption index 

of upto a < However, we note that using a more finely tuned setting of the constant e in the proof of 


Theorem 15 and a more careful proof using tight tail inequalities for chi-squared distributions [7], we can 


achieve a better corruption level tolerance of a < A. 


□ 


D Proof of Theorem 


Theorem 5. Let X = [xi,... ,x„] S be the given data matrix and y = A^w* -|- b be the corrupted 

output with ||b||o < a ■ n. Let X satisfy the SSC and SSS properties at level 7 with constants and A.^ 
respectively (see Definition [^. Let Algorithmbe executed on this data with the GD update (Algorithm]^ 
with the thresholding parameter set to 13 > a and the step length set to 77 = —. If the data satisfies 

max |? 7 y^A^, 1 — 77 Ai_^} < |, then after t = O ^log iterations. Algorithm 

solution w* i.e. Ilw* — w*|L < e. 


obtains an e-accurate 


Proof. Let r* = y — A w* be the vector of residuals at time t and Ct = Xs^Xg^. We have 


w 


‘+1 = w‘ + 77 • As,4^ = w* -h ?7 • As,(yst - 


The thresholding step ensures that 


J+i 


— Il’'st^ll 2 Claim 
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and use (3 > a) which implies 


J+i 


< 


.t+1 

‘^MD, 


where CRt+i = St+i\St are the corrupted recoveries and MD^+i = SA\St+i are the clean points missed out 
from detection. Note that |CRt+i| < a ■ n and |MDt_|_i| < f3 ■ n. Since bs, = 0 and MDj+i C 5*, we get 


bcRt+i + ^CRt+i (w* - w + ) 


< 


A 


MD. 


n(w*-w + ) 


Using the SSS conditions and the fact that ||bst_,.i II 2 = ||bst+i\s. II 2 gives us 


||bst+i II2 = ||bcRt+i II2 < (\/^ + V^) 

Now, using the expression for gives us 


d+i| 


< 2443 I 


,t+i| 


w — w 




< ||(/-77Ct)(w* - w*)|| +77||As,bs, 


We will bound the two terms on the right hand separately. We can bound the second term easily as 

»?l|2fstbsj|2 < ?7v^l|bstll2 < I|bstll2: 

since Ijb ^,||q < a • n. For the first term we observe that for 77 < we have 

\\I-r]Ct\\ 2 = sup |l - ?7 • v^Civ] = sup {l - 77 • v^Ctv} < 1 - 77Ai_/3, 
vGSp-1 vGSp-1 


which we can use to bound 


7 + 11 


< (1-77Ai_^) 


l|bsj 
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This gives us, for r] = 

||bst+i||2 - 2\/^| 


,t+i| 


<21- 






(P) 


Ai_/5 

(Q) 


|bst 


For Gaussian designs and small enough /?, we can show (Q) < ^ as we did in Theorem H To bound (P), we 
use the lower bound on Xi-p given by Theorem 15 and use the following tighter upper pound for Ai_^: 


Ai _/3 < ^(1 - f3) + 36^6/3(1 - f3) log n + O (^Jnp + n log 


The above bound is obtained similarly to the one in Theorem 15 but uses the identity (^) = < 

( \ n—k 

j for values of fc > n/2 instead. For small enough /3 and n = VL («;^(S)(p + log j)), we can then 
show (P) < j as well. Let dtj := i/u ||w* — w‘|j2 + || 5 st II • Using elementary manipulations and the fact that 
\/A^ > n {\/n), we can then show that 

^i+i < 3/4 • 4-4. 

Thus, int = O ^log ^^|lw*||2 + iterations of the algorithm, we arrive at an e-optimal solution i.e. 

||w* — w *||2 < e. A similar argument holds true for sub-Gaussian designs as well. □ 


E Proof of Theorem [6] 

Theorem 6. Suppose Algorithm is executed on data that allows Algorithms and a convergence rate of 
77 fc and rycD respectively. Suppose we have 2 • rjFc ‘ Vgb < 1- Then for any interleaviMs of the FG and GD 
steps that the policy may enforce, after t = O ^log iterations. Algorithm 4 ensures an e-optimal 

solution i.e. ||w‘ — w*|| < e. 

Proof. Our proof shall essentially show that the FC and GD steps do not undo the progress made by the 
other if executed in succession and if 2 • rypc • Vgd < 1; actually ensure non-trivial progress. Let 

= l|bsJl2 

= V^||w* - w*|| -h Ilbsj|2 

denote the potential functions used in the analyses of the FG and GD algorithms before. Then we will show 
below that if the FC and GD algorithms are executed in steps t and t -|- 1 then we have 

< 2 • 77fc • T?GD • 

Alternatively, if the GD and FC algorithms are executed in steps t and t -|- 1 respectively, then 

'L™ < 2 • r]Fc ■ ilGD ■ 

Thus, if algorithm executes the FC step at the time step t, then it would at least ensure < (2 • 77 fc ■ Vgb Y^^' 
(similarly if the last step is a GD step). Since both the FC and GD algorithms ensure ||w* — w *||2 < e 

for t > O ^log ^be claim would follow. 

We now prove the two claimed results regarding the two types of interleaving below 
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1. FC —^ GD 

The FC step guarantees ||bsj_^J |2 < Tyre • ||bsj| as well as ||w‘+^ “'''^*112 — , whereas the 

GD step guarantees < rycD • Together these guarantee 

y/n ||w*+^ - W *||2 + ||bst+ 2||2 - ■ V^||w*+^ - W*\\^ + ||bSt+i ||2 

< 2 • rjFc ■ VGI} ■ I!bsjl2 

Since ^/n ||w*+^ ~ '''^*||2 — yields the result. 

2. GD—^FC 

The GD step guarantees < rycD • whereas the FC step guarantees < ?7 fc • ||bst_,.i || 

as well as ||w‘+^ ~ '''^*||2 — ■ Together these guarantee 

Vn||w*+2 - w^ll^ ^ 2?7 fc ||bst+i||2 

< 2 • ryFC • ??GD ■ 

where the second step follows from the GD step guarantee since i/n ||w*+^ — w *||2 > 0. 

This finishes the proof. □ 


F Proof of Theorem [9] 


Theorem 9. Let X = [xi,... ,x„] S be the given data matrix and y = X'^w* + b be the corrupted 

output with ||w*|jg < s* and ||b||o < a ■ n. Let Algorithmbe executed on this data with the IHT update 
from [5] and thresholding parameter set to (3 > a. Let Eg be an invertible matrix such that Eg satisfies 
the SRSC and SRSS properties at level (7,2s + s*) with constants a(g, 2 s+s*) and L(g., 2 s+s*) respectively 

(see Definition 8 ) for s > 32 with 7 = 1 — /3. If A also satisfies < 1 then after 


t = o 




iterations 


, Algorithm 


obtains an e-accurate solution w‘ i.e. ||w* — w*|j 2 < e. In 


particular, if A is sampled from a Gaussian distribution A/’(0, E) and n > D ((2s + s*) logp + log |), then 
for all values of a < /3 < we can guarantee recovery as H-w* — w*II 2 < e. 

Proof. We first begin with the guarantee provided by existing sparse recovery techniques. The results of 
[5], for example, indicate that if the input to the algorithm indeed satisfies the RSC and RSS properties 
at the level (1 — ft, 2s + s*) with constants a 2 s+s* and L 2 s+s* for s > 32 ^, then in time r = 

^ (bfew ' algorithm [5J Algorithm 1 ] outputs an updated model w*+^ that satisfies 


|w‘+^||g < s, as well as 


T „,*+! I 


ys, - Ag^w 


< ||yst - X]w* 


P- 


We will set p later. Since the SRSC and SRSS properties ensure the above and y = A ' w* + b, this gives us 


Ag^(w‘"'"^ — w *)||2 < 2 (w*“'"^ — w*)' Ag^bgj + p = 2 (w*“'"^ — w*)^A^ 


Stns.^Stns. 


■ P, 


since bs = 0 for any set S' n S'* = (/). We now analyze the two sides separately below using the SRSC and 

~ _1 /o 

SRSS properties below. For any S C [n], denote A 5 := Eg ' A. 


|Aj(w‘+i - w*)||J = AjEy"(w‘+i - w*) 


^ Q(l_y3,s+s*) 


Eg (w — W ) 


A 


St ns. 


(w‘+^ - w*) = 


A 


St ns. "2^0 




(w*+-' - w*) 


< 


■>(0,8 + 3*) 


G/2/,„t+l 


(w‘+^ - W*) 
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Now, if ||w‘+i 


— w 


I2 ^ 


then 


- w*) > ^/Xnuni^o) ■ e- This give us 


y\l/2 

^0 


(w‘+^ 


^ 2\/T(/3,S + S*) IIj^ _ II , P 

„ - -iPStns. II2 + 7- 

2 ^(1—^,s+s’“) ^(1-/3,s+s"*) 

_ 2-v/T(/3,s+s*) II P 

— - -I|bsjl2 H-/ ■ ' -• 

^(1~/5,s+s*) e • \/Amin(^oj ‘ ^(1-/9,s+s*) 


We note that although we declared the SRSC and SRSS properties for the action of matrices on sparse 
vectors (such as w* — w‘+^), we instead applied them above to the action of matrices on sparse vectors 
transformed by — w*+^)). Since need not be sparse even if v is sparse, this appears 

to pose a problem. However, all we need to resolve this is to notice that the proof technique of Theorem [T8| 
which would be used to establish the SRSC and SRSS properties, holds in general for not just the action of a 
matrix on the set of sparse vectors, but on vectors in the union of any fixed set of low dimensional subspaces. 

More specifically, we can modify the RSC and RSS properties (and by extension, the SRSC and SRSS 
properties), to requiring that the matrix X act as an approximate isometry on the following set of vectors 

:= |v : V = Eg for some v' S We refer the reader to the work of [2] which describes this 

technique in great detail. Proceeding with the proof, the assurance of the thresholding step, as used in the 
proof of Theorem along with a straightforward application of the (modified) SRSS property gives us 


|bst+i|L < 


^CRt+i(w*+^ - W*) 


+ 


- W*) 


-^CRt+i^O 


^1/2 .■ t-\-\ * 

En (W — W 


+ 


X 


T 

MDt+i^O 


U/2 / t+l 


(w*“^ — W*) 


< 2JL 


^(/3,s+s*) 




(w‘+^ - W*) 


^(1-/3,s+s*) 


‘^P\/^{P,S+S*) 


min (Eo) • a(i_^_s+s*) 


Thus, whenever ||w‘+^ — w* > e, in successive steps, |lb 5 j |2 undergoes a linear decrease. Denoting 

T] := ^^1£i£±£E_ get 


|bs,+J2 ^ P' 




s+s*) 


l-rf 

^ ^ / C • yj Aniin(Eo) * 


and using 


Ey^(w*-w*) > \/Amin(Eo) ||w‘ - w *||2 gives us 


|w*+^ — w’ll^ < - , _ 

V -^min (Eg) • «( 

1—/3,s+s*) 


2'\/'^(/3,s+s*) 


|bst+i II + 


P 




2\/T(/3, 


s+s*) 


Aniin(Eo) ■ 0^1 — ^ 

mp 


\A min (Eo) • o;(i_;3^s+s*) 


■^min (Eg) • a(i_^_s+s 


where we have assumed that < 9/10, something that we shall establish below. Note that 

Amin (Eg) > 0 since E is assumed to be invertible. In the random design settings we shall consider, we 
also have ■\/^(g,.+3 > - _ q /-EV Then setting p < ,Tg 2 , s+s*) proves the conver- 

■^A„i„(So)'a(i_,3,a + a*) VV"/ y r, J 

gence result. 

As before, we can use the above result to establish sparse recovery guarantees in the statistical setting for 
Gaussian and sub-Gaussian design models. If our data matrix X is generated from a Gaussian distribution 
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Af{0, S) for some invertible S, then the results in Theorem 18 can be used to establish that S satisfies 

the SRSC and SRSS properties at the required levels and that for a < and n> Q ((2s + s*) logp + log |), 
we have ri = ^ 9/10. 

Thus, the above result can be applied with Sq = 51 to get convergence guarantees in the general Gaussian 
setting. We note that the above analysis can tolerate the same level of corruption as Theorem]^ and thus, 
we can improve the noise tolerance level to a < ^ here as well. We also note that these results can be 
readily extended to the sub-Gaussian setting as well. □ 


G Robust Statistical Estimation 

This section elaborates on how results on the convergence guarantees of our algorithms can be used to give 
guarantees for robust statistical estimation problems. We begin with a few definition of sampling models 
that would be used in our results. 

Definition 12. A random variable a: S K is called sub-Gaussian if the following quantity is finite 

sup 

P>1 

Moreover, the smallest upper bound on this quantity is referred to as the sub-Gaussian norm of x and denoted 
as ||a;||^^. 

Definition 13. A vector-valued random variable x € is called sub-Gaussian if its unidimensional 
marginals (x, v) are sub-Gaussian for all v G 5'^“^. Moreover, its sub-Gaussian norm is defined as fol¬ 
lows 

11^1^2 ■= II(^>V)||^^ 

vGSp-1 

We will begin with the analysis of Gaussian designs and then extend our analysis for the class of general 
sub-Gaussian designs. 

Lemma 14. Let X G be a matrix whose columns are sampled i.i.d from a standard Gaussian distri¬ 

bution i.e. Xi ~ JV{0,I). Then for any e > 0, with probability at least 1 — 5, X satisfies 

Sn,a.^{XX^) < n -b (1 - 2e)~^ ^cnp -\- c'n log 

Smin{XX^) > n - (1 - 2e)~^ ^cnp -\- c'n log 
where c = 24e^ log | and d = 24e^ . 

Proof. We will first use the fact that X is sampled from a standard Gaussian to show that its covariance 
concentrates around identity. Thus, we first show that with high probability, 

-nJ||2 < ei 

for some ei < 1. Doing so will automatically establish the following result 

n-ei< s„,iniXX^) < Sma^XX^) < n -b ei. 

Let A := XX^ — I. We will use the technique of covering numbers m to establish the above. Let 
CP~^(e) C be an e cover for i.e. for all u G there exists at least one v G such 

that ||u —v ||2 < e. Standard constructions [CT see Lemma 5.2] guarantee such a cover of size at most 
(l -b 1)^ < (f )^- Now for any u G and v G such that ||u — v ||2 < e, we have 

— v^ 24 v| < |u^ 24 (u — v)| -b |v^ 24 (u — v)| < 2e ||v4||2, 
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which gives us 


\XX^-nl\\ <{l-2e)-'^ ■ sup 

vGCP-i(e) 




Now for a fixed v € 5”" the random variable ||X ' vjl^ is distributed as a x^(n) distribution with n degrees 
of freedom. Using Lemma 20 we get, for any /i < 1, 


F 


> fj,n 


< 2 exp — min 


. f /rn 

— mm < — 


i 24ne2 ’ 4^36 


< 2 exp 


24e2 ) ■ 


Setting fj,"^ = c - ^ + c' ■ '"^here c = 24e^ log | and c' = 24e^, and taking a union bound over all ^(e), 
we get 

\ 


sup 

veCP-i(e) 


F 


> cnp + dn log - 


<2|5V»p|- 


24e2 ) - 


< 5. 


This implies that with probability at least 1 — <5, 


IF^ — n/jl^ < (1 — 2 e) Y cnp + c'nlog -, 
which gives us the claimed bounds on the singular values of XX^. 


□ 


Theorem 15. Let X € be a matrix whose columns are sampled i.i.d from a standard Gaussian 

distribution i.e. ~ A/'(0, 1). Then for any 7 > 0, with probability at least 1 — 5, the matrix X satisfies the 
SSC and SSS properties with constants 


A, 


Gauss 


< jn [ 1 + 3e. / 6 log - + O \ np + n\og 


7 


1 


L Gauss 


> n — (1 — 7 )n ( 1 + 3e. /h log 


1-7 


- U I y np + n log - I . 


Proof. For any fixed S € S^, Lemma 14 guarantees the following bound 


Smi,AXsXg) < 7 n+(l- 2 e) ^^cynp + Uynlog 
Taking a union bound over S.y and noting that (^) < for all 1 < /c < n, gives us 


e 2 

A..y < yn + (1 — 2e)~^\lc'ynp + log —h Uynlog - 

70 


< 7 n 1 + (1 — 2 e) ./c'log - + (1 — 2 e) W cynp + c'yn log- 


which finishes the first bound after setting e = 1/6. For the second bound, we use the equality 

XsXj = XXT - XgXj, 
which provides the following bound for X.y 


> S„,in(XXT) - sup XtX^ = S^iniXX^) - Ai_ 


TgSi-.~, 


7- 


Using Lemma to bound the first quantity and the first part of this theorem to bound the second quantity 
gives us, with probability at least 1 — 5, 


Xj > n — • 


7 'n + (1 - 2 e) ^ ^d log - (1 - 2 e) ^ + vF) c-np + dn log ?, 

where 7 ' = 1 — 7 . This proves the second bound after setting e = 1/6. 


□ 
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We now extend our analysis to the class of isotropic subGaussian distributions. We note that this analysis 
is without loss of generality since for non-isotropic sub-Gaussian distributions, we can simply use the fact 
that Theorem [S] can admit whitened data for calculation of the SSC and SSS constants as we did for the 
case of non-isotropic Gaussian distributions. 

Lemma 16. Let X G be a matrix with columns sampled from some sub-Gaussian distribution with 

sub-Gaussian norm K and covariance E. Then, for any 5 > 0, with probability at least 1 — 5, each of the 
following statements holds true: 

^max(^^ ) ^ '5niax(^) ‘ ^ “t” “f ty^ 

SminiXX^) > Ami„(E) ■ n - Ck ' 


where t = y ^log|, and ck,Ck are absolute constants that depend only on the sub-Gaussian norm K of 
the distribution. 


Proof. Since the singular values of a matrix are unchanged upon transposition, we shall prove the above 
statements for X^. The beneht of this is that we get to work with a matrix with independent rows, so that 
standard results can be applied. The proof technique used in PH Theorem 5.39] (see also Remark 5.40 (1) 
therein) can be used to establish the following result: with probability at least 1 — 5, with t set as mentioned 
in the theorem statement, we have 


WX' -E 


< Ck 




This implies that for any v € S'p we have 


— ||X^v||^ — v^Ev 

71 '' ^ 


-v^XX^v - v^Ev 
n 


< 


-XX ' V - Ev 
n 


<CksI- 
n 


The results then follow from elementary manipulations and the fact that the singular values and eigenvalues 
of real symmetric matrices coincide. □ 

Theorem 17. Let X S IRp^" be a matrix with columns sampled from some sub-Gaussian distribution with 
sub-Gaussian norm K and covariance E. Let ck, Ck and t be fixed to values as required in Lemma \T^ Note 
that ck and Ck are absolute constants depend only on the sub-Gaussian norm K of the distribution. Let 
7 G (0,1] be some fixed eonstant. Then, with we have the following: 

^ subGaussiK ,^) < ^ log -U + Ck- ^ + ty/fl. 


Furthermore, fix any e G (0,1) and let ^ be a value in (0,1) satisfying the following 


7 > 1 — min 


e ■ Amin(S) 

Aniax(E) 


exp 


(^1 + IR-i 


^2 \2 
' '^min 


(s) 


where W-i{-) is the lower branch of the real valued restriction of the Lambert W function. Then we have, 
with the same confidence, 

yubGaussiK,^) > _ 2g) . .n-CK(l+ x/T^) - 2tyCi 

Proof. The first result follows from an application of Lemma pH a union bound over sets in as well as the 
bound (2) < for all 1 < fc < n which puts a bound on the number of sparse sets as log |5.^| < 7 -nlog C 
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For the second result, we observe that — XgX^, so that Sniin(-^s-^s) — — 

Smax(-^S^J)- This gives us 

iuf Sniin(^S^5 ) — ^raini^^ ) SUp Sniax(^S^5)- 
■ 565 ^ SeSi-^ 

Using Lemma and the first part of this result gives us 

^n| Smi„(XsXj) > Amin(S) • 71 - Cif • ^/]m - 

Amax(S)(l - 7 ) + \ -—- log . ^ ) n - CkV{^- 7)p^ - tVn 

\ ck 1 - 7 / 

= (^A„,in(S) - A^ax(S)(l - 7 ) - y^log n 

-Ck (1 + \/l - 7 ) - 2tv^ 

> (1 - 2e) • Amin(S) -n-CK (^ + \/l - 7 ) “ 2t^/n, 

where the last step follows from the assumptions on 7 and by noticing that it suffices to show the following 
two inequalities to establish the last step 

1 - Ai„ax(S){l - 7 ) < e • Amin(U) 

2 . (1 - 7 ) log < cxe^ ■ A^i„(E) 

The first part gives us the condition 7 > 1 — in a straightforward manner. For the second part, 

^max 1“^ j 

denote v = ■ A^ijj(E). Note that for u > 1, all values of 7 S (0,1] satisfy the inequality. 

Otherwise we require the use of the Lambert W function (also known as the product logarithm function). 
This function ensures that its value W{z) for any z > —1/e satisfies z = W{z)e^^^\ In our case, making a 
change of variable (1 — 7 ) = gives us the inequality (77— > —vje. Note that since u < 1 in this case, 

—vje G (—1/e, 0) i.e. a valid value for the Lambert W function. However, (—1/e, 0) is also the region in 
which the Lambert W function is multi-valued. Taking the worse bound for 7 by choosing the lower branch 
lF_i(-) gives us the second condition 7 > 1 — exp ^1 -|- lF_i 'Nninl^l ^y □ 

It is important to note that for any —1/e < z < 0, we have exp (1 -|- lF_i(z)) > 0 which means that the 
bounds imposed on 7 by Theorem always allow a non-zero fraction of the data points to be corrupted 
in an adversarial manner. However, the exact value of that fraction depends, in a complicated manner, on 
the sub-Gaussian norm of the underlying distribution, as well as the condition number and the smallest 
eigenvalue of the second moment of the underlying distribution. 

We also note that due to the generic nature of the previous analysis, which can handle the entire class of 
sub-Gaussian distributions, the bounds are not as explicitly stated in terms of universal constants as they 
are for the standard Gaussian design setting (Theorem [T^ . 

We now establish that for a wide family of random designs, the SRSG and SRSS properties are satisfied 
with high probability as well. For sake of simplicity, we will present our analysis for the standard Gaussian 
design. However, the results would readily extend to general Gaussian and sub-Gaussian designs using 
techniques similar to Theorem [l7l 

Theorem 18. Let X G he a matrix whose columns are sampled i.i.d from a standard Gaussian 

distribution i.e. ~ Then for any 7 > 0 and s < p, with probability at least 1 — <5, the matrix X 

satisfies the SRSC and SRSS properties with constants 



T Gauss 


< 777 


(^l-b3e 



-f O 


I ns + n log - 
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^ Gauss 
“(7.S) 


> n — (1 — 7 )n ( 1 + 3ew61og 


1-7 


- n 


I ns + n log - 


Proof. The proof of this theorem proceeds similarly to that of Theorem |15| Hence, we simply point out the 
main differences. First, we shall establish, that for any e > 0, with probability at least 1 — 5, X satisfies the 
RSC and RSS properties at level s with the following constants 


Lg < n + (1 — 2e) \ bns + b'nlog - 


a., > n — (1 — 2e) ^ \l bns + b'nlog-, 


where b = 24e^ log ^ and b' = 24e^. To do so we notice that the only change needed to be made would be 
in the application of the covering number argument. Instead of applying the union bound over an e-cover 
of SP~^, we would only have to consider an e-cover of the set of all s-sparse unit vectors in 
p-dimensions. A straightforward calculation shows us that 




-11 


< 


e 


< 


Sep 

es 


Thus, setting fj? = b ■ ^b' ■ where b = 24e^ log ^ and b' = 24e^, we get 


sup 

vGcr 


||Xv ||2 -n 


> \ bns + b'n log 


< A 


which establishes the required RSC and RSS constants for X. Now, moving on to the SRSS constant, it 
follows simply by applying a union bound over all sets in S.y much like in Theorem 15 One can then proceed 
to bound the SRSC constant in a similar manner. 

We note that the nature of the SRSC and SRSS bounds indicate that our Torrent-FC algorithm in 
the high dimensional sparse recovery setting has noise tolerance properties, characterized by the largest 
corruption index a that can be tolerated, identical to its low dimnensional counterpart - something that 
Theorem states explicitly. □ 


H Supplementary Results 

Claim 19. Given any vector v S M", let a G Sn be defined as the permutation that orders elements of v 
in descending order of their magnitudes i.e. |^^cr(i)| > |'^o'( 2 )| > > |i'(T(n)|- Abr any 0 < p < g < 1, let 

Si G Sq be an arbitrary set of size q • n and S 2 = {o'(f) :n — p-n + l<i< n}. Then we have Hvs^H^ < 
f||vsj|^<||v5j|^ 

Proof. Let S 3 = {cr{i) : n — q-n + l<i<n} and S '4 = :n — q-n + l<i<n — p - n}. Clearly, we 

have IIV 53 II 2 — IKS 1 II 2 since S 3 contains the smallest q ■ n elements (by magnitude). Now we have IIVS 3 II 2 = 
ll''^S 2 ll 2 + ll''^S' 4 ll 2 - Moreover, since each element of S 4 is larger in magnitude than every element of S 2 , we 
have 

This gives us 

IIVS2II2 = IIVS3II2 - IIVS4II2 < IIVS3II2 - ||^ lIvSallL 

which upon simple manipulations, gives us the claimed result. □ 
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Lemma 20. Let Z be distributed according to the chi-squared distribution with k degrees of freedom i.e. 
Z ~ Then for all t>0, 


\Z — k\ > t] < 2 exp ( — 


( t 

{ 24fce2 ’ 4^36 


Proof. This lemma requires a proof structure that traces several basic results in concentration inequalities for 
sub-exponential variables [H Lemma 5.5, 5.15, Proposition 5.17]. The purpose of performing this exercise 
is to explicate the constants involved so that a crisp bound can be provided on the corruption index that 
our algorithm can tolerate in the standard Gaussian design case. 

We first begin by establishing the sub-exponential norm of a chi-squared random variable with a single 
degree of freedom. Let X ~ Then using standard results on the moments of the standard normal 

distribution gives us, for all p > 2, 

(E|Xr)Vp = ((2p - 1 )!!)Vp = 

Thus, the sub-exponential norm of X is upper bounded by •\/3/2. By applying the triangle inequality, 
we obtain, as a corollary, an upper bound on the sub-exponential norm of the centered random variable 

Y = X-le.s\\Y\\^^<2\\X\\^^<VS._ 

Now we bound the moment generating function of the random variable Y. Noting that ET = 0, we have, 
for any |A| < 


Eexp(AF) 


^ ^ ^ E(Ay)^ 

^ n \ 


< 1 q. ^ 


<1-1- y^(-\/3e|A|)'^ < 1 + 6e^A^ < exp(6e^A^). 

q=2 


Note that the second step uses the sub-exponentially of Y, the third step uses the fact that q\ > {q/eY, and 
the fourth step uses the bound on |A|. Now let Xi, X 2 ,... Xk be k independent random variables distributed 
as X^(l)- Then we have Z ^ Using the exponential Markov’s inequality, and the independence of 

the random variables Xi gives us 


¥[Z - k>t]=F ]jEexp(A(Xi - 1)). 

For any |A| < the above bounds on the moment generating function give us 


k 

¥[Z — k>t] < JJexp(6e^A^) = exp(—At -I- Gfce^A^). 

i=l 


Choosing A = min < 


\ 2V3: 


12ke^ 


|, we get 


¥[Z-k>t] 


< exp 


min 


_ 

24fce2 ’ 4v^e J / 


Repeating this argument gives us the same bound for ¥\k — Z >t]. This completes the proof. 


□ 
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I Supplementary Experimental Results 


n = 2000 alpha = 0.25 sigma = 0.2 

Ir 


0.5 


■■■TORRENT-FC 
♦TORRENT-HYB 
, —L1-DALM 



1000 


p = 500 n = 2000 alpha = 0.25 


*TORRENT-FC 

♦TORRENT-HYB 



Sigma {White Noise) 


(b) 


p = 500 alpha = 0.25 sigma = 0.2 P = 500 n = 2000 alpha = 0.25 sigma = 0.2 




(d) 


Figure 3: (a), (b), (c) Variation of recovery error with varying p,a and n. Torrent was found to outperform 
DALM-Li in all these settings, (d) Recovery error as a function of runtime for various state-of-the-art Li solvers as 
indicated in m- 
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